SFC: A Macro-mapping

Antoine Godin (Kingston University)

AB-SFC Workshop - University of Leeds - May 19

Aims: SFC modelling and national accounts

Network of authors in 2013

Network of authors in 2013, source: Caverzasi and Godin (2015)

Network of authors in 2013, source: Caverzasi and Godin (2015)

Network of authors in 2016

Network of authors in 2016, source: author’s computation

Network of authors in 2016, source: author’s computation

National Accounts and Stock-Flow Consistent Modelling

Stock-Flow Accounting

Balance Sheets

Non-financial corporations

Assets Liabilities
Produced non-financial asset 1857180 0
Non-produced non financial assets 0 0
Currency and deposits 546337 0
Securities other than shares 65083 384631
Loans 262326 965308
Shares and other equity 837888 2475658
Insurance technical reserves 4029 1056253
Other accounts receivable/payable 29976 50912
Net Worth 0 -1386059

Financial corporations

Assets Liabilities
Produced non-financial asset 113216 0
Non-produced non financial assets 0 0
Currency and deposits 3129884 5141459
Securities other than shares 3520842 1829957
Loans 3418831 1553526
Shares and other equity 2986514 2046996
Insurance technical reserves 1296602 3844915
Other accounts receivable/payable 6069641 6031487
Net Worth 0 15088

Goverment

Assets Liabilities
Produced non-financial asset 1857180 0
Non-produced non financial assets 0 0
Currency and deposits 546337 0
Securities other than shares 65083 384631
Loans 262326 965308
Shares and other equity 837888 2475658
Insurance technical reserves 4029 1056253
Other accounts receivable/payable 29976 50912
Net Worth 0 -1386059

Balance Sheets in SFC

Example

HHs Firms Gov. Banks C. B. Sum
Capital +Kh +Kf +K
Money +Hh +Hb -H 0
Bills +Bh -Bs +Bb +Bcb 0
Loans -Lh -Lf +L 0 0 0
Equities +Ef -Ef 0 0 0 0
Equities +Eb 0 -Eb 0 0 0
Net worth -NWh -NWf -NWb -NWg 0 -K
Sum 0 0 0 0 0 0

Sectorial accounts (from Eurostat)

Example for households in the UK in 2014

# Selecting the flows
names<-c("B5G","D5","D61","D62","D7","D8","B6G","P3","B8G","P5G","D9","NP","B9")

# Obtaining the data
EZdata_raw = pdfetch_EUROSTAT("nasa_10_nf_tr", UNIT="CP_MNAC",NA_ITEM=names, GEO="EU28",
                              SECTOR=c("S14_S15"), TIME="2014")

# Transforming the data into a data.frame
EZdata<-as.data.frame(EZdata_raw)

# Automatic procedure to remove the non-interesting bit of the colnames
coln<-colnames(EZdata)
newcoln<-c()
HHdata<-c()

for(i in 1:length(coln)){
  name<-coln[i]
  tname<-strsplit(name,"\\.")[[1]]
  newname<-paste(tname[3:4],collapse=".")
  # If the column contains only NA, remove it from the dataset
  if(!is.na(EZdata[16,i])){
    newcoln<-c(newcoln,newname) 
    HHdata<-c(HHdata,EZdata[16,i])
  }
}

# Creating a new dataset with only values 2014
HHdata<-as.data.frame(t(HHdata))
colnames(HHdata)<-newcoln

# Creating the aggregates
HHdata_1<-as.data.frame(c(HHdata$PAID.B5G,-HHdata$PAID.D5,-HHdata$PAID.D61+HHdata$RECV.D61,
                          +HHdata$RECV.D62-HHdata$PAID.D62,-HHdata$PAID.D7+HHdata$RECV.D7,paste(HHdata$PAID.B6G,"]"),-HHdata$PAID.P3,-HHdata$PAID.D8+HHdata$RECV.D8,paste(HHdata$PAID.B8G,"]"),-HHdata$PAID.P5G,-HHdata$PAID.D9+HHdata$RECV.D9,
                          -HHdata$PAID.NP,HHdata$PAID.B9))
colnames(HHdata_1)<-"Households"
rownames(HHdata_1)<-c("Total Income","Taxes","Social Contributions","Social Benefits",
                      "Other transfers","[Gross Disposable Income","Consumption","Adjustments in Pensions","[Gross Savings","Gross Capital Formation","Capital Transfer",
                      "Net Non-Produced NF Assets","Net Lending Position")
kable(HHdata_1)
Households
Total Income 9881408
Taxes -1455366
Social Contributions -2422712
Social Benefits 2637790
Other transfers 116967
[Gross Disposable Income 8758087 ]
Consumption -8008156
Adjustments in Pensions 200560
[Gross Savings 950492 ]
Gross Capital Formation -710480
Capital Transfer 23045
Net Non-Produced NF Assets 6961
Net Lending Position 270017

Transaction flow matrix

  1. The transactions flow matrix consists of three separate parts.
    • On the top rows of the matrix, you will have output expressed as expenditures, which by definition is given by \(Y = C + I + G (+ X - M)\)
    • The matrix should clearly identify the consumption and investment by the sectors in your model.
  2. The second part of the transactions flow matrix outlines output using an income approach.
    • Depending on how disaggregated your model is and how many assets you have included in your model, this part may include various sources of income for your sectors

Assume a closed economy,

Example of Transaction Flow Matrix

Transaction Flow Matrix, source: G&L 2007

Transaction Flow Matrix, source: G&L 2007

Transaction flow matrix, part 2

Income ? Expenditures ? Change in assets/liabilities =0

Full integration matrix

Example of Full Integration Matrix

Full Integration Matrix, source: G&L 2007

Full Integration Matrix, source: G&L 2007

Data

Eurostat

Data Search in Eurostat

Data Search in Eurostat

Data Table in Eurostat

Data Table in Eurostat

Data Details

Data Details

Looking at sectoral Accounts

url: http://ec.europa.eu/eurostat/web/sector-accounts

Sectoral Accounts webpage

Sectoral Accounts webpage

Sectoral Accounts Countries Data

Sectoral Accounts Countries Data

pdfetch: getting data automatically

Example 1: Net lending per sector, UK

# Specifying the name of the flows of interests
names<-c("B9")
# Downloading the data by specifying the various filters
UKdata_raw = pdfetch_EUROSTAT(flowRef = "nasa_10_nf_tr", 
        UNIT="CP_MNAC",NA_ITEM=names, GEO="UK", DIRECT="PAID",
        SECTOR=c("S11","S12","S13","S14_S15","S2"))
# Transforming the obtained data into a data.frame
UKdata<-as.data.frame(UKdata_raw)
# Setting readable names
colnames(UKdata)<-c("NFC","FC","Govt","HH","RoW")

# Matricial plot
matplot(as.numeric(substr(row.names(UKdata),1,4)),UKdata,lwd=2,type="l",lty=1,
        ylab="",xlab="",col=2:6)
# Adding the horizontal line
abline(h=0,col=1)
# Adding a grid
grid()
# Adding a legend
legend("bottomleft",col=c(2:6),lwd=2,lty=1,legend=c("NFC",
        "FC","Govt","HH","RoW"),bty='n')

Example 2: Household income statement from 2014

# Selecting the flows
names<-c("B5G","D5","D61","D62","D7","D8","B6G","P3","B8G"
            ,"P5G","D9","NP","B9")
# Obtaining the data
EZdata_raw = pdfetch_EUROSTAT("nasa_10_nf_tr", 
            UNIT="CP_MNAC",NA_ITEM=names, GEO="EU28",
            SECTOR=c("S14_S15"), TIME="2014")
# Transforming the data into a data.frame
EZdata<-as.data.frame(EZdata_raw)

# Automatic procedure to remove the non-interesting bit
# of the colnames
coln<-colnames(EZdata)
newcoln<-c()
HHdata<-c()
for(i in 1:length(coln)){
  name<-coln[i]
  tname<-strsplit(name,"\\.")[[1]]
  newname<-paste(tname[3:4],collapse=".")
# If the column contains only NA, remove it from dataset
  if(!is.na(EZdata[16,i])){
    newcoln<-c(newcoln,newname) 
    HHdata<-c(HHdata,EZdata[16,i])
  }
}

# Creating a new dataset with only values 2014
HHdata<-as.data.frame(t(HHdata))
colnames(HHdata)<-newcoln
# Creating the aggregates
HHdata_1<-as.data.frame(c(HHdata$PAID.B5G,-HHdata$PAID.D5,
            -HHdata$PAID.D61+HHdata$RECV.D61,
            +HHdata$RECV.D62-HHdata$PAID.D62,
            -HHdata$PAID.D7+HHdata$RECV.D7,HHdata$PAID.B6G))
colnames(HHdata_1)<-"Households"
rownames(HHdata_1)<-c("Total Income","Taxes"
        ,"Social Contributions","Social Benefits",
        "Other transfers","Gross Disposable Income")
kable(HHdata_1)
Households
Total Income 9881408
Taxes -1455366
Social Contributions -2422712
Social Benefits 2637790
Other transfers 116967
Gross Disposable Income 8758087
HHdata_2<-as.data.frame(c(HHdata$PAID.B6G,-HHdata$PAID.P3,
        -HHdata$PAID.D8+HHdata$RECV.D8,HHdata$PAID.B8G))
colnames(HHdata_2)<-"2014"
rownames(HHdata_2)<-c("Gross Disposable Income",
    "Consumption","Adjustments in Pensions","Gross Savings")
kable(HHdata_2)
2014
Gross Disposable Income 8758087
Consumption -8008156
Adjustments in Pensions 200560
Gross Savings 950492
HHdata_3<-as.data.frame(c(HHdata$PAID.B8G,-HHdata$PAID.P5G
        ,-HHdata$PAID.D9+HHdata$RECV.D9,-HHdata$PAID.NP,
        HHdata$PAID.B9))

colnames(HHdata_3)<-"2014"
rownames(HHdata_3)<-c("Gross Savings",
        "Gross Capital Formation","Capital Transfer",
        "Net Non-Produced NF Assets","Net Lending Position")
kable(HHdata_3)
2014
Gross Savings 950492
Gross Capital Formation -710480
Capital Transfer 23045
Net Non-Produced NF Assets 6961
Net Lending Position 270017

PKSFC package

Installing the dependent libraries and the package

You need to install all the required libraries This is for traditional libraries

install.packages("expm")
install.packages("igraph")

For non-conventional libraries, such as the one need to visualize Direct Acyclical Graphs (DAG), you need to do the following

source("http://bioconductor.org/biocLite.R")
biocLite("Rgraphviz")

Finally you can then download the PKSFC package from github and install it locally

install.packages("path/PKSFC_1.5.tar.gz", repos = NULL, 
                                 type="source")

Testing

Now we’re ready to load the package:

library(PKSFC)
## Loading required package: expm
## Loading required package: Matrix
## 
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
## 
##     expm
## Loading required package: igraph
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
  1. Load SIM (download SIM.txt from Github)
sim<-sfc.model("Models/SIM.txt",modelName="SIMplest model 
        from chapter 3 of Godley and Lavoie (2007)")
  1. Simulate the model
datasim<-simulate(sim)
  1. replicate figure 3.2 of page 73.
matplot(sim$time,datasim$baseline[,c("Yd","C")],type="l", 
    xlab="", ylab="",lwd=2,lty=c(2,3))
lines(sim$time,vector(length=length(sim$time))
    +datasim$baseline["2010","C"], lwd=2)
legend(x=1970,y=50,lwd=2, legend=c("Disposable Income",
    "Consumption","Steady State"),lty=c(2,3,1), bty="n")

How does it work?

Source code of SIM

#1. EQUATIONS
C = C
G = G
T = T
N = N
Yd = W*N - T
T = theta*W*N
C = alpha1*Yd + alpha2*H_h(-1)
H = H(-1) + G - T
H_h = H_h(-1) + Yd - C
Y = C + G
N = Y/W
#2. PARAMETERS
alpha1=0.6
alpha2=0.4
theta=0.2
#EXOGENOUS
G=20
W=1
#INITIAL VALUES
H=0
H_h=0
#3. Timeline
timeline 1945 2010

A few important points regarding the model source code:

Internal representation

print(sim)
## $name
## [1] "SIMplest model \n\t\tfrom chapter 3 of Godley and Lavoie (2007)"
## 
## $simulated
## [1] FALSE
## 
## $variables
##       name     initial value description    
##  [1,] "Yd"     NA            "1. EQUATIONS" 
##  [2,] "T"      NA            ""             
##  [3,] "C"      NA            ""             
##  [4,] "H"      "0"           ""             
##  [5,] "H_h"    "0"           ""             
##  [6,] "Y"      NA            ""             
##  [7,] "N"      NA            ""             
##  [8,] "alpha1" "0.6"         "2. PARAMETERS"
##  [9,] "alpha2" "0.4"         ""             
## [10,] "theta"  "0.2"         ""             
## [11,] "G"      "20"          "EXOGENOUS"    
## [12,] "W"      "1"           ""             
## 
## $endogenous
##      name  initial value lag description   
## [1,] "Yd"  NA            "0" "1. EQUATIONS"
## [2,] "T"   NA            "0" ""            
## [3,] "C"   NA            "0" ""            
## [4,] "H"   "0"           "1" ""            
## [5,] "H_h" "0"           "1" ""            
## [6,] "Y"   NA            "0" ""            
## [7,] "N"   NA            "0" ""            
## 
## $equations
##      endogenous value equation                 description   
## [1,] "Yd"             "W*N-T"                  "1. EQUATIONS"
## [2,] "T"              "theta*W*N"              ""            
## [3,] "C"              "alpha1*Yd+alpha2*H_h_1" ""            
## [4,] "H"              "H_1+G-T"                ""            
## [5,] "H_h"            "H_h_1+Yd-C"             ""            
## [6,] "Y"              "C+G"                    ""            
## [7,] "N"              "Y/W"                    ""            
## 
## $time
##  [1] 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958
## [15] 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972
## [29] 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986
## [43] 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000
## [57] 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010
## 
## $matrix
##     Yd T C H H_h Y N
## Yd   0 1 0 0   0 0 1
## T    0 0 0 0   0 0 1
## C    1 0 0 0   0 0 0
## H    0 1 0 0   0 0 0
## H_h  1 0 1 0   0 0 0
## Y    0 0 1 0   0 0 0
## N    0 0 0 0   0 1 0
## 
## $blocks
## $blocks[[1]]
## [1] 2 3 4 6 7 1 5
## 
## 
## attr(,"class")
## [1] "sfc"

Output data structure

The Gauss Seidel Algorithm

  1. Select initial values \(x^0\)
  2. While \(k<maxIter\) & \(\delta < tolValue\)
  1. For each \(i=1,...,n\): \[x_i^{k+1}=\frac{1}{a_{ii}}\left( b_i-\sum^{i-1}_{j=1}a_{ij} x_j^{k+1}-\sum_{j=i+1}^n a_{ij}x_j^k\right)\]
  2. Compute \(\delta\): \[\delta = \frac{x^{k+1}-x^k}{x^k}\]

System of (in)dependent equations

See Fenell et. al (2016)

Direct Acyclic graphs

simex<-sfc.model("Models/SIMEX.txt",modelName="SIMplest model with expectation")
layout(matrix(c(1,2),1,2))
plot.dag(sim,main="SIM" )
plot.dag(simex,main="SIMEX" )

In the case of a more complex model - Chapter 6

ch6 <- sfc.model("Models/ch6.txt",modelName="Chapter6_openmodel")
graphs = generate.DAG.collapse( adjacency = ch6$matrix )
par(mfrow = c(1,3))
# first plot the orgianl grpah
plot_graph_hierarchy( graphs$orginal_graph, main = "orginal graph" )
# plot hte nodes that form the strongly connected compoent
plot_graph_hierarchy( graphs$SCC_graph, main = "SCC nodes" )
# plot the result ing DAG when we take the condensation of the graph
plot_graph_hierarchy( graphs$DAG, main = "resulting DAG" )

Systems of dependent vs independent equations

#Doing the simulations
datasimex<-simulate(simex)
init = datasimex$baseline[66,]
simex<-sfc.addScenario(simex,"G",25,1960,2010,init)
datasimex<-simulate(simex)
datasim<-simulate(sim)
init = datasim$baseline[66,]
sim<-sfc.addScenario(sim,"G",25,1960,2010,init)
datasim<-simulate(sim)

#Plotting it all
layout(matrix(c(1,2),1,2))
matplot(sim$time,datasim$scenario_1[,c("H","C","Yd")],
    lty=1:3,type="l",xlab="",ylab="",main="SIM")
legend(x=1944,y=130,legend=c("Wealth","Consumption",
    "Disposable Income"), lty=1:3,bty="n")
matplot(simex$time,datasimex$scenario_1[,c("H","C","Yd",
    "Yd_e")],lty=1:4,type="l",xlab="",ylab="",main="SIMEX")
legend(x=1944,y=130,legend=c("Wealth","Consumption","Disposable Income",
   "Expecetd Disposable Income"),lty=1:4,bty="n")

Computational implications

Let’s see how much time it takes to run sim:

ptm <- proc.time()
data1<-simulate(sim)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],
                        "seconds"))
## [1] "Elapsed time is  3.35 seconds"

Now lets play with some of the parameters of the simulate function:

  1. tolValue
ptm <- proc.time()
data2<-simulate(sim,tolValue = 1e-3)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],"seconds"))
## [1] "Elapsed time is  0.190000000000001 seconds"
  1. maxIter
ptm <- proc.time()
data3<-simulate(sim, maxIter=10)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],"seconds"))
## [1] "Elapsed time is  0.210000000000001 seconds"

Observing the results of the three simulations

kable(round(t(cbind(data1$baseline[c(1,2,20,40,66),c("Y")],
        data2$baseline[c(1,2,20,40,66),c("Y")],
        data3$baseline[c(1,2,20,40,66),c("Y")])),digits=3))
1945 1946 1964 1984 2010
NA 38.462 96.957 99.892 99.999
NA 38.464 96.874 99.577 99.911
NA 34.006 94.965 99.672 99.991

Block Gauss-Seidel

The order of equations matters, if first compute variables that do not depend on current period, this speeds the process. Define blocks of equation independent from the others.

print(simex$blocks)
## [[1]]
## [1] 7
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 6
## 
## [[4]]
## [1] 8
## 
## [[5]]
## [1] 2
## 
## [[6]]
## [1] 1 4
## 
## [[7]]
## [1] 5

Simulation of SIMEX

ptm <- proc.time()
dataex<-simulate(simex,tolValue = 1e-10)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3]
                        ,"seconds"))
## [1] "Elapsed time is  0.16 seconds"

Results for SIMEX

kable(round(t(dataex$baseline[c(1,2,20,40,66),c("G","Y",
    "T","Yd","Yd_e","C","H","H_h")]),digits=3))
1945 1946 1964 1984 2010
G 20 20 20.000 20.000 20
Y NA 20 98.559 99.983 100
T NA 4 19.712 19.997 20
Yd 0 16 78.847 79.987 80
Yd_e NA 0 78.559 79.983 80
C NA 0 78.559 79.983 80
H 0 16 78.847 79.987 80
H_h 0 16 78.847 79.987 80

Output data structure

kable(dataex$baseline)
Yd T C H H_h Y Yd_e N alpha1 alpha2 theta G W iter block 1 iter block 2 iter block 3 iter block 4 iter block 5 iter block 6 iter block 7
1945 0.00000 NA NA 0.00000 0.00000 NA NA NA 0.6 0.4 0.2 20 1 0 0 0 0 0 0 0
1946 16.00000 4.00000 0.00000 16.00000 16.00000 20.00000 0.00000 20.00000 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1947 28.80000 7.20000 16.00000 28.80000 28.80000 36.00000 16.00000 36.00000 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1948 39.04000 9.76000 28.80000 39.04000 39.04000 48.80000 28.80000 48.80000 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1949 47.23200 11.80800 39.04000 47.23200 47.23200 59.04000 39.04000 59.04000 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1950 53.78560 13.44640 47.23200 53.78560 53.78560 67.23200 47.23200 67.23200 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1951 59.02848 14.75712 53.78560 59.02848 59.02848 73.78560 53.78560 73.78560 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1952 63.22278 15.80570 59.02848 63.22278 63.22278 79.02848 59.02848 79.02848 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1953 66.57823 16.64456 63.22278 66.57823 66.57823 83.22278 63.22278 83.22278 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1954 69.26258 17.31565 66.57823 69.26258 69.26258 86.57823 66.57823 86.57823 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1955 71.41007 17.85252 69.26258 71.41007 71.41007 89.26258 69.26258 89.26258 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1956 73.12805 18.28201 71.41007 73.12805 73.12805 91.41007 71.41007 91.41007 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1957 74.50244 18.62561 73.12805 74.50244 74.50244 93.12805 73.12805 93.12805 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1958 75.60195 18.90049 74.50244 75.60195 75.60195 94.50244 74.50244 94.50244 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1959 76.48156 19.12039 75.60195 76.48156 76.48156 95.60195 75.60195 95.60195 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1960 77.18525 19.29631 76.48156 77.18525 77.18525 96.48156 76.48156 96.48156 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1961 77.74820 19.43705 77.18525 77.74820 77.74820 97.18525 77.18525 97.18525 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1962 78.19856 19.54964 77.74820 78.19856 78.19856 97.74820 77.74820 97.74820 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1963 78.55885 19.63971 78.19856 78.55885 78.55885 98.19856 78.19856 98.19856 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1964 78.84708 19.71177 78.55885 78.84708 78.84708 98.55885 78.55885 98.55885 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1965 79.07766 19.76942 78.84708 79.07766 79.07766 98.84708 78.84708 98.84708 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1966 79.26213 19.81553 79.07766 79.26213 79.26213 99.07766 79.07766 99.07766 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1967 79.40970 19.85243 79.26213 79.40970 79.40970 99.26213 79.26213 99.26213 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1968 79.52776 19.88194 79.40970 79.52776 79.52776 99.40970 79.40970 99.40970 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1969 79.62221 19.90555 79.52776 79.62221 79.62221 99.52776 79.52776 99.52776 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1970 79.69777 19.92444 79.62221 79.69777 79.69777 99.62221 79.62221 99.62221 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1971 79.75821 19.93955 79.69777 79.75821 79.75821 99.69777 79.69777 99.69777 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1972 79.80657 19.95164 79.75821 79.80657 79.80657 99.75821 79.75821 99.75821 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1973 79.84526 19.96131 79.80657 79.84526 79.84526 99.80657 79.80657 99.80657 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1974 79.87621 19.96905 79.84526 79.87621 79.87621 99.84526 79.84526 99.84526 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1975 79.90096 19.97524 79.87621 79.90096 79.90096 99.87621 79.87621 99.87621 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1976 79.92077 19.98019 79.90096 79.92077 79.92077 99.90096 79.90096 99.90096 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1977 79.93662 19.98415 79.92077 79.93662 79.93662 99.92077 79.92077 99.92077 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1978 79.94929 19.98732 79.93662 79.94929 79.94929 99.93662 79.93662 99.93662 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1979 79.95944 19.98986 79.94929 79.95944 79.95944 99.94929 79.94929 99.94929 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1980 79.96755 19.99189 79.95944 79.96755 79.96755 99.95944 79.95944 99.95944 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1981 79.97404 19.99351 79.96755 79.97404 79.97404 99.96755 79.96755 99.96755 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1982 79.97923 19.99481 79.97404 79.97923 79.97923 99.97404 79.97404 99.97404 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1983 79.98338 19.99585 79.97923 79.98338 79.98338 99.97923 79.97923 99.97923 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1984 79.98671 19.99668 79.98338 79.98671 79.98671 99.98338 79.98338 99.98338 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1985 79.98937 19.99734 79.98671 79.98937 79.98937 99.98671 79.98671 99.98671 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1986 79.99149 19.99787 79.98937 79.99149 79.99149 99.98937 79.98937 99.98937 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1987 79.99319 19.99830 79.99149 79.99319 79.99319 99.99149 79.99149 99.99149 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1988 79.99456 19.99864 79.99319 79.99456 79.99456 99.99319 79.99319 99.99319 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1989 79.99564 19.99891 79.99456 79.99564 79.99564 99.99456 79.99456 99.99456 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1990 79.99652 19.99913 79.99564 79.99652 79.99652 99.99564 79.99564 99.99564 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1991 79.99721 19.99930 79.99652 79.99721 79.99721 99.99652 79.99652 99.99652 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1992 79.99777 19.99944 79.99721 79.99777 79.99777 99.99721 79.99721 99.99721 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1993 79.99822 19.99955 79.99777 79.99822 79.99822 99.99777 79.99777 99.99777 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1994 79.99857 19.99964 79.99822 79.99857 79.99857 99.99822 79.99822 99.99822 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1995 79.99886 19.99971 79.99857 79.99886 79.99886 99.99857 79.99857 99.99857 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1996 79.99909 19.99977 79.99886 79.99909 79.99909 99.99886 79.99886 99.99886 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1997 79.99927 19.99982 79.99909 79.99927 79.99927 99.99909 79.99909 99.99909 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1998 79.99942 19.99985 79.99927 79.99942 79.99942 99.99927 79.99927 99.99927 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1999 79.99953 19.99988 79.99942 79.99953 79.99953 99.99942 79.99942 99.99942 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2000 79.99963 19.99991 79.99953 79.99963 79.99963 99.99953 79.99953 99.99953 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2001 79.99970 19.99993 79.99963 79.99970 79.99970 99.99963 79.99963 99.99963 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2002 79.99976 19.99994 79.99970 79.99976 79.99976 99.99970 79.99970 99.99970 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2003 79.99981 19.99995 79.99976 79.99981 79.99981 99.99976 79.99976 99.99976 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2004 79.99985 19.99996 79.99981 79.99985 79.99985 99.99981 79.99981 99.99981 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2005 79.99988 19.99997 79.99985 79.99988 79.99988 99.99985 79.99985 99.99985 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2006 79.99990 19.99998 79.99988 79.99990 79.99990 99.99988 79.99988 99.99988 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2007 79.99992 19.99998 79.99990 79.99992 79.99992 99.99990 79.99990 99.99990 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2008 79.99994 19.99998 79.99992 79.99994 79.99994 99.99992 79.99992 99.99992 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2009 79.99995 19.99999 79.99994 79.99995 79.99995 99.99994 79.99994 99.99994 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2010 79.99996 19.99999 79.99995 79.99996 79.99996 99.99995 79.99995 99.99995 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2

Checking the number of iteractions for SIM

kable(round(t(cbind(
    data1$baseline[c(1,2,20,40,66),c("iter block 1")],
    data2$baseline[c(1,2,20,40,66),c("iter block 1")],
    data3$baseline[c(1,2,20,40,66),c("iter block 1")]
    )),digits=3))
1945 1946 1964 1984 2010
0 293 218 175 119
0 89 15 2 1
0 10 10 10 10

Checking the number of iteractions for simex, no simulate parameters

kable(round(t(dataex$baseline[c(1,2,20,40,66),
    c("iter block 1","iter block 2","iter block 3",
    "iter block 4","iter block 5","iter block 6",
    "iter block 7")]),digits=3))
1945 1946 1964 1984 2010
iter block 1 0 2 2 2 2
iter block 2 0 2 2 2 2
iter block 3 0 2 2 2 2
iter block 4 0 2 2 2 2
iter block 5 0 2 2 2 2
iter block 6 0 2 2 2 2
iter block 7 0 2 2 2 2

Visualisation

Visualisation - Sankey Diagram

Flows in the UK Economy 2014, source: Burgess et al. (2016)

Flows in the UK Economy 2014, source: Burgess et al. (2016)

Visualisation - Graph representation

Example with models SIM-SIMEX

Graph representation of SIM and SIMEX, source: author’s computation

Graph representation of SIM and SIMEX, source: author’s computation

Shock propagation in more complex models

Graph representation of Burgess et al. (2016)

Graph representation of Burgess et al. (2016)

Shock to CAR - Direct

Graph representation of Burgess et al. (2016)

Graph representation of Burgess et al. (2016)

Shock to CAR - Lag 1: Interest rates

Graph representation of Burgess et al. (2016)

Graph representation of Burgess et al. (2016)

Shock to CAR - Lag 2: Income and profits

Graph representation of Burgess et al. (2016)

Graph representation of Burgess et al. (2016)

Dynamical interaction with models

antoinegodin.shinyapps.io/SIMEX

Shiny application

Shiny application

Model with Portfolio choice

Balance Sheet

Households Production Government Central Bank Sum
Money +H -H 0
Bills +Bh -B +Bcb 0
Net worth -V +V 0
Sum 0 0 0 0 0

Transaction Flow Matrix

Households Production Government Central Bank Sum
Consumption -C +C 0
Gov. Exp. +G -G 0
Income = GDP +Y -Y 0
Interests +r(-1)*Bh(-1) -r(1)*B(-1) +r(-1)*Bcb(-1) 0
CB profits +r(-1)*Bcb(-1) -r(1)*Bcb(-1) 0
Taxes -T +T 0
Change in Money -\(\Delta\) H +\(\Delta\) H 0
Change in Bills -\(\Delta\) Bh +\(\Delta\) B -\(\Delta\) Bcb 0
Sum 0 0 0 0 0

Behavioural Equations (from PC.txt file)

#######Determination of output - eq. 4.1
y = cons + g
#######Disposable income - eq. 4.2
yd = y - t + r(-1)*b_h(-1)
#######Tax payments - eq. 4.3
t = theta*(y + r(-1)*b_h(-1))
#######Wealth accumulation - eq. 4.4
v = v(-1) + (yd - cons)
#######Consumption function - eq. 4.5
cons = alpha1*yd + alpha2*v(-1)
#######Cash money - eq. 4.6
h_h = v - b_h  
#######Demand for government bills - eq. 4.7
b_h = v*(lambda0 + lambda1*r - lambda2*(yd/v))
#######Supply of government bills - eq. 4.8
b_s = b_s(-1) + (g + r(-1)*b_s(-1)) - (t + r(-1)*b_cb(-1))
#######Supply of cash - eq. 4.9
h_s = h_s(-1) + b_cb - b_cb(-1)
#######Government bills held by the central bank - eq. 4.10
b_cb = b_s - b_h
#######Interest rate as policy instrument - eq. 4.11
r = r_bar

Portfolio Equations

How to allocate wealth between bills and money?

  1. Quantity theory of money: links money balances to the flow of income
  2. Liquidity preferences: money balances some proportion of total wealth

Merging the two into Tobin’s portfolio equations:

\[\begin{align} \frac{H_h}{V}&=(1-\lambda_0)-\lambda_1\cdot r+\lambda_2\cdot\left(\frac{YD}{V}\right)\tag{4.6A}\\ \frac{B_h}{V}&=\lambda_0+\lambda_1\cdot r-\lambda_2\cdot\left(\frac{YD}{V}\right)\tag{4.7}\\ H_h&=V-B_h\tag{4.6} \end{align}\]

Setting up the environment

Before doing any modelling, we need to load the package in the R environment.

library(PKSFC)

Then, you need to download the attached ‘PC.txt’ file and save it in the folder of your choice. Make sure to set the working directory where you saved the downloaded file. In command line this looks like this but if you use Rstudio, you can use the graphical interface as well (Session>Set Working Directory>Choose Directory)

setwd("pathToYourDirectory")

Loading the model

The first thing to do is to load the model and check for completeness.

pc<-sfc.model("Models/PC.txt",modelName="Portfolio Choice Model")
pc<-sfc.check(pc,fill=FALSE)

Let’s have a look at the graphical representation of model PC. In order to do so, we need to load the Rgraphviz library:

library("Rgraphviz")

We can now look at the graph of model PC:

plot.dag(pc,main="PC" )

You can see that there is a cycle in the graph, implying that GDP, taxes, disposable income and consumption are determined all together (and that they fully adapt to any shock applied to the economy). While this is a mathematical property of the system of equations representing the economy we wish to model, it has an economic meaning and you want to be sure that this is what you believe to be the best representation of what you have in mind.

We are now ready to simulate the model

datapc<-simulate(pc)

Expectations and random shocks

The first of experiment is meant to observe the buffer role of money, in case of random shocks applied to expected disposable income. To do these, we need to modify slightly model pc and change the equations determining consumption, demand for bonds and money, and the expectations on income and wealth. We will call the new model pcRand.

pcRand<-sfc.editEqus(pc,list(list(var="cons",eq="alpha1*yde+alpha2*v(-1)"),
      list(var="b_h",eq="ve*(lambda0 + lambda1*r - lambda2*(yde/ve))")))
pcRand<-sfc.addEqus(pcRand,list(
  list(var="yde",eq="yd(-1)*(1+rnorm(1,sd=0.1))",
         desc="Expected disposable income depending on random shocks"),
  list(var="h_d",eq="ve-b_h",desc="Money demand"),
  list(var="ve",eq="v(-1)+yde-cons",desc="Expected disposable income")))
pcRand<-sfc.editVar(pcRand,var="yd",init=86.48648)
pcRand<-sfc.check(pcRand)

You probably have noticed that the yde equation contains the R function rnorm which allows you to extract a random number from a normal distribution centred around 0, with standard variation of 0.1. The package indeed allows you to use any R function such as min, max or rnorm. However, because of the way the Gauss-Seidel algorithm works, we are facing a problem. Indeed, as the rnorm function extract a number in each iteration, this will mean that the algorithm cannot converge since the difference between each iteration of the algorithm depends of the difference between each extraction. This is why we need to restrict the number of iteration of the Gauss-Seidel.

Before simulating the model, let’s have a look at how the graph of the model has changed:

plot.dag(pcRand,main="PC Random" )

You can now see that the cycle observed in the original model has disappeared. Indeed, consumption doesn’t depend on disposable income any more but on expected disposable income. Let’s now simulate the model.

datapcRand<-simulate(pcRand,maxIter=2)

This replicates figure 4.1, page 110

plot(pcRand$time,datapcRand$baseline[,"h_h"],type="l",xlab="",ylab="",lty=1,lwd=2,
     ylim=c(1*min(datapcRand$baseline[,c("h_h","h_d")],na.rm=T),
            1.2*max(datapcRand$baseline[,c("h_h","h_d")],na.rm=T)))
lines(pcRand$time,datapcRand$baseline[,"h_d"],lty=2,lwd=2)
legend(x=1950,y=1.2*max(datapcRand$baseline[,c("h_h","h_d")],na.rm=T),
       legend=c("Money held","Money demand"),lty=c(2,1),lwd=2,bty="n")

This graph highlights the buffer role of certain stocks in PK-SFC models. Indeed, because expectations are incorrect or because demand might not be equal to supply in any market, at least one stock will not be equal to the targeted level. As highlighted by Foley (1975), in a model without perfect foresight you need a buffer stock in order to obtain equilibrium between demand and supply. The role of buffer stocks in PK-SFC model is thus fundamental and is at the hart of the approach used by Godley (1999) in his seven unsustainable processes. It is by observing the dynamics of certain stock-flow norms that you are able to observe the unsustainable processes evolving in an economy, because stocks precisely absorb disequilibrium.

Interest rates impacts

The graphs presented on the paper are based on the model PCEX, which includes expectation on disposable income and wealth. We will first create the model before simulating it.The shock represents an 100 basis point increase in interest rates.

pcex<-sfc.editEqus(pc,list(
  list(var="cons",eq="alpha1*yde+alpha2*v(-1)"),
  list(var="b_h",eq="ve*(lambda0 + lambda1*r - lambda2*(yde/ve))")))
pcex<-sfc.addEqus(pcex,list(
  list(var="yde",eq="yd(-1)",desc="Expected disposable income"),
  list(var="h_d",eq="ve-b_h",desc="Money demand"),
  list(var="ve",eq="v(-1)+yde-cons",desc="Expected disposable income")))
pcex<-sfc.editVar(pcex,var="yd",init=86.48648)
pcex<-sfc.check(pcex)
init = datapc$baseline[56,]
pcex<-sfc.addScenario(model=pcex,vars="r_bar",values=0.035,inits=1960,
                                            ends=2010,starts=init)
datapcex<-simulate(pcex)

This replicates plot figure 4.3. p. 112

time2=c(1950:2000)
plot(time2, datapcex$scenario[as.character(time2),"h_h"]/
            datapcex$scenario[as.character(time2),"v"],
         type="l",xlab="",ylab="",lty=1,lwd=2)
par(new=T)
plot(time2,datapcex$scenario[as.character(time2),"b_h"]/
            datapcex$scenario[as.character(time2),"v"],
     lty=2, lwd=2,type="l",axes=F,ylab="",xlab="")
axis(4,pretty(c(0.750, 1.1*max(datapcex$scenario[,"b_h"]/
            datapcex$scenario[,"v"],na.rm=T))))
legend(x=1970,y=0.78,legend=c("Share of Bills",
        "Share of Money"),lty=c(2,1),lwd=2,bty="n")
grid()

No surprise here, the change in return rates coming from the shock led to a re-allocation between Bills and Money, in favour of the former.

This replicates fig 4.4 . p 113

time2=c(1950:2000)
matplot(time2,datapcex$scenario[as.character(time2),c("yd","cons")]
                ,type="l",xlab="",ylab="",lty=1:2,col=1,lwd=2)
legend(x=1970,y=88,legend=c("Disposable Income","Consumption"),
             lwd=c(2,2),lty=c(1,2),bty="n")
grid()

These results are more surprising as increased interest rates lead to an increase both in disposable income and consumption. To understand this, one has to bear in mind that increase in interest rates leads to an increase in government spending and thus to an increase in the fiscal stance which determines the steady state:

\(Y^\star=\frac{G+r \cdot B^{\star}_h \cdot (1-\theta)}{\theta}\)

Endogenous propensities to consume

We change the propensities to consume to incorporate the fact that they might be impacted by interest rates and thus have

\[\begin{equation} \alpha_1=\alpha_10-\iota\cdot r_{-1} \tag{4.32} \end{equation}\]

This piece of code shows how to implement this

pc_end<-sfc.addEqus(pcex,list(
  list(var="alpha1",eq="alpha10-iota*r(-1)")))
pc_end<-sfc.addVars(pc_end,list(
  list(var="alpha10",init="0.7",
         desc="Endoegnous propensity to consume - autonomous term"),
  list(var="iota",init="10",
         desc="Endoegnous propensity to consume - interest rate impact")))
pc_end$scenarios<-NULL
pc_end<-sfc.check(pc_end)
datapc_end<-simulate(pc_end)
init = datapc_end$baseline[56,]
pc_end<-sfc.addScenario(model=pc_end,vars="r_bar",values=0.035,
                                                inits=1960,ends=2010,starts=init)
datapc_end<-simulate(pc_end)

This replicates plot figure 4.9. p. 123

time2=c(1950:2000)
matplot(time2,
     datapc_end$scenario_1[as.character(time2),c("v","yd","cons","y")],
     type="l",xlab="",ylab="",lty=1:4,lwd=2,col=1)
grid()
legend(x=1950,y=100,legend=c("Households wealth","Disposable income",
            "Consumption","GDP"),lty=1:4,lwd=2,bty="n")

This shows that now the short-term dynamics display the Keynesian paradox of thrift but the long-run steady state still shows a positive impact of interests rates on GDP. The short-run recession is governed by the decrease in the propensity to consumed and is slowly compensated by the increase in wealth (and the increased public debt as a counterpart) leading to the new steady state. The steady state is still driven by the fiscal stance.

This replicates plot figure 4.10. p. 123

timelag=c(1950:2000)
time2=c(1951:2001)
plot(time2,
     datapc_end$scenario_1[as.character(time2),c("g")]+
            datapc_end$scenario_1[as.character(timelag),c("r")]*
            datapc_end$scenario_1[as.character(timelag),c("b_s")]-
            datapc_end$scenario_1[as.character(timelag),c("r")]*
            datapc_end$scenario_1[as.character(timelag),c("b_cb")],
     type="l",xlab="",ylab="",lty=1:4,lwd=2,col=1,ylim=c(20.8,23))
lines(time2,datapc_end$scenario_1[as.character(time2),"t"],type="l",xlab="",
            ylab="",lty=2,lwd=2)
grid()
legend("topleft",legend=c("Government expenditures plus net debt services",
            "Tax revenues"),lty=1:4,lwd=2,bty="n")

Debt to GDP at the steady state

\[\begin{equation} \frac{V^\star}{Y^\star}=\frac{\frac{1-\alpha_1}{\alpha_2}}{1+\left[\frac{\theta}{1-\theta}\right]-r\cdot \left[(\lambda_0+\lambda_1\cdot r)\cdot \frac{1-\alpha_1}{\alpha_2}-\lambda_2\right]}\tag{4.33} \end{equation}\]

The Maastricht treaty: The reference values referred to […] are: - 3% for the ratio of the planned or actual government deficit to gross domestic product at market - 60% for the ratio of government debt to gross domestic product at market prices.

Adding the Maastricht treaty: 1 Rational Government

#Simulate it
datapcex<-simulate(pcex)
#Computing the targeted tax rate
init <- as.data.frame(t(datapcex$baseline[66,]))
alpha1<-init$alpha1
alpha2<-init$alpha2
alpha<-(1-alpha1)/alpha2
r<-init$r
lambda0<-init$lambda0
lambda1<-init$lambda1
lambda2<-init$lambda2
lambda<-lambda0+lambda1*r
M<-0.6
newTheta <- (alpha/M+r*(lambda*alpha-lambda2)-1)/(alpha/M+r*(lambda*alpha-lambda2))

Adding a scenario

Changing the value of the tax rate to its optimal value

pcex<-sfc.addScenario(pcex,c("theta"),c(newTheta),c(1960),c(2010),init)
#Simulation
datapcex<-simulate(pcex)

Doing some plots

  1. debt to GDP
#Plots the results
plot(pcex$time,datapcex$scenario_2[,"v"]/datapcex$scenario_2[,"y"],
    type="l",xlab="",ylab="",xlim=range(pcex$time),main="Debt to GDP")

  1. GDP, consumption and disposable income
plot(pcex$time,datapcex$scenario_1[,"y"],type="l",xlab="",ylab="",
    xlim=range(pcex$time),ylim=range(datapcex$scenario_1[,c("cons","yd","y")],na.rm=T))
lines(pcex$time,datapcex$scenario_1[,"cons"],lty=2)
lines(pcex$time,datapcex$scenario_1[,"yd"],lty=3)
legend(x=1970,y=110,legend=c("GDP","Consumption","Disposable Income"),lty=c(1,2,3),bty="n")

A (somewhat) Rational Government

Updating the model to a growth government expenditure model

pcexgr<-sfc.addEqu(pcex,"g","g(-1)*(1+grg)","government expenditures")
pcexgr<-sfc.editVar(pcexgr,var="theta",init=newTheta)
pcexgr<-sfc.addVar(pcexgr,var="grg",init=0,desc="Government expenditure growth")
pcexgr$scenarios<-NULL
pcexgr<-sfc.addScenario(pcexgr,list(c("grg")),list(c(0.03)),1960,1982,init)
pcexgr<-sfc.check(pcexgr)


#Simulation
datapcex2<-simulate(pcexgr)

Plots the results

  1. Deficit to GDP
plot(pcexgr$time[2:66],(datapcex2$scenario_1[2:66,"v"]-datapcex2$scenario_1[1:65,"v"])
        /datapcex2$scenario_1[2:66,"y"],type="l",xlab="",ylab="",
        xlim=range(pcexgr$time),main="Deficit to GDP")

  1. GDP, consumption and disposable income
plot(pcexgr$time,datapcex2$scenario_1[,"y"],type="l",xlab="",ylab="",
    xlim=range(pcexgr$time),ylim=range(datapcex2$scenario_1[,c("cons","yd","y")],na.rm=T))
lines(pcexgr$time,datapcex2$scenario_1[,"cons"],lty=2)
lines(pcexgr$time,datapcex2$scenario_1[,"yd"],lty=3)
legend(x=1965,y=110,legend=c("GDP","Consumption","Disposable Income"),lty=c(1,2,3),bty="n")

A (more) Rational Government

Updating the model to endogenous tax rate model

theta<-0.2
pcexgr2<-sfc.addEqu(pcex,"g","g(-1)*(1+grg)","government expenditures")
pcexgr2<-sfc.addEqu(pcexgr2,"theta","theta(-1)+dtheta","Tax rate")
pcexgr2<-sfc.addVar(pcexgr2,var="grg",init=0,desc="Government expenditure growth")
pcexgr2<-sfc.addVar(pcexgr2,var="dtheta",init=0,desc="Tax rate variation")
pcexgr2$scenarios<-NULL
pcexgr2<-sfc.addScenario(pcexgr2,list(c("grg","dtheta")),list(c(0.03,(newTheta-theta)/24))
                    ,1960,1982,init)
pcexgr2<-sfc.check(pcexgr2)

#Simulation
datapcex3<-simulate(pcexgr2)

Plots the results

  1. Deficit to GDP
plot(pcexgr$time[2:66],(datapcex3$scenario_1[2:66,"v"]-datapcex3$scenario_1[1:65,"v"])
        /datapcex3$scenario_1[2:66,"y"],type="l",xlab="",ylab="",xlim=range(pcexgr$time),
        main="Deficit to GDP")

  1. GDP, consumption and disposable income
plot(pcexgr$time,datapcex3$scenario_1[,"y"],type="l",xlab="",ylab="",
    xlim=range(pcexgr$time),ylim=range(datapcex3$scenario_1[,c("cons","yd","y")],na.rm=T))
lines(pcexgr$time,datapcex3$scenario_1[,"cons"],lty=2)
lines(pcexgr$time,datapcex3$scenario_1[,"yd"],lty=3)
legend(x=1965,y=110,legend=c("GDP","Consumption","Disposable Income"),lty=c(1,2,3),bty="n")

Liquidity Preferences and model LP

Value of perpetuity, interest rates, exepcted returns and capital gains

The value of a financial asset is supposed to reflect the net present value of future cash flows. Thus for a long-term bond who is never redeemed, the price \(p_{bL}\) of the bond is given by

\(p_{bL} = \sum \frac{1}{(1+r_{bL})^t} = \frac{1}{r_{bL}}\)

Return rate of an asset is equal to both the yield (interests, dividends) and the capital gains normalized to the nominal value of the asset. In the case of a long-term bond:

\(Rr_{bL} = r_{bL}(-1)+\frac{\Delta p_{bL}}{p_{bL}(-1)} = \frac{1+p_{bL}-p_{bL}(-1)}{p_{bL}(-1)}\)

However, what matter for the households when making their decision is the expected price of bonds \(p_{bL}^e\), given the current price of bonds. This leads to the pure expected rate of return.

\(PERr_{bL} = r_{bL} + \frac{p^e_{bL}-p_{bL}}{p_{bL}}\)

Because expectation are not followed by everyone, households might incorporate a weight \(\xi\) into their expectation formation. This leads to the expected return rate

\(ERr_{bL} = r_{bL} + \xi \frac{p^e_{bL}-p_{bL}}{p_{bL}}\)

The change in prices of bonds leads to capital gains which is one of two parts of the change in nominal value of financial assets. More precisely:

\[\begin{align*} \Delta (p_{BL}\cdot BL) &=\left(p_{bL} \cdot BL\right) - \left(p_{bL}(-1) \cdot BL(-1)\right)\\ &= p_{BL} \cdot BL-p_{BL}\cdot BL(-1) + p_{BL}\cdot BL(-1)-p_{bL}(-1) \cdot BL(-1)\\ &= p_{BL} \left( BL-BL(-1)\right) + BL(-1) \left(p_{BL}-p_{BL}(-1)\right)\\ &= p_{BL} \cdot \Delta BL+ BL(-1)\cdot \Delta p_{BL} \end{align*}\]

Thus total changes in nominal values of an asset is equal to the quantity of assets bought (or sold) at the current price (\(p_{BL} \cdot \Delta BL\)) and to the capital gains (or losses) made on the stock of assets held at the beginning of the period (\(BL(-1)\cdot \Delta p_{BL}\))

Ostergaard diagram

Ostergaard diagram

Balance Sheet

Households Production Government Central Bank Sum
Money +H -H 0
Bills +Bh -B +Bcb 0
Bonds +BL.pbL -BL.pbL 0
Net worth -V +V 0
Sum 0 0 0 0 0

Transaction Flow Matrix

Households Production Government Central Bank Sum
Consumption -C +C 0
Gov. Exp. +G -G 0
Income = GDP +Y -Y 0
Interests on bills +rb(-1)*Bh(-1) -rb(1)*B(-1) +rb(-1)*Bcb(-1) 0
Interests on bonds +BL(-1) -BL(-1) 0
CB profits +r(-1)*Bcb(-1) -r(1)*Bcb(-1) 0
Taxes -T +T 0
Change in Money -\(\Delta\) H +\(\Delta\) H 0
Change in Bills -\(\Delta\) Bh +\(\Delta\) B -\(\Delta\) Bcb 0
Change in Bonds -\(\Delta\) BL.pbL +\(\Delta\) BL.pbL 0
Sum 0 0 0 0 0
Memo: Capital gains -\(\Delta\) pbL.BL(-1) +\(\Delta\) pbL.BL(-1) 0

Equation list

# MODEL
# Determination of output - eq. 5.1
y = cons + g
# Regular disposable income - eq. 5.2
yd_r = y - t + r_b(-1)*b_h(-1) + bl_h(-1)
# Tax payments - eq. 5.3
t = theta*(y + r_b(-1)*b_h(-1) + bl_h(-1))
# Wealth accumulation - eq. 5.4
v = v(-1) + (yd_r - cons) + cg
# Capital gains on bonds - eq. 5.5
cg = (p_bl-p_bl(-1))*bl_h(-1)
# Consumption function - eq. 5.6
cons = alpha1*yd_r_e + alpha2*v(-1)
# Expected wealth - eq. 5.7
v_e = v(-1) + (yd_r_e - cons) + cg
# Cash money - eq. 5.8
h_h = v - b_h - p_bl*bl_h
# Demand for cash - eq. 5.9
h_d = v_e - b_d - p_bl*bl_d
# Demand for government bills - eq. 5.10
b_d = v_e*(lambda20 + lambda22*r_b - lambda23*er_rbl - lambda24*(yd_r_e/v_e))
# Demand for government bonds - eq. 5.11
bl_d = v_e*(lambda30 - lambda32*r_b + lambda33*er_rbl - lambda34*(yd_r_e/v_e))/p_bl
# Bills held by households - eq. 5.12
b_h = b_d
# Bonds held by households - eq. 5.13
bl_h = bl_d
# Supply of government bills - eq. 5.14
b_s = b_s(-1) + (g + r_b(-1)*b_s(-1) + bl_s(-1)) - (t + r_b(-1)*b_cb(-1)) -p_bl*(bl_s-bl_s(-1))
# Supply of cash - eq. 5.15
h_s = h_s(-1) + b_cb - b_cb(-1)
# Government bills held by the central bank - eq. 5.16
b_cb = b_s - b_h
# Supply of government bonds - eq. 5.17
bl_s = bl_h
# Expected rate of return on bonds - eq. 5.18
er_rbl = r_bl+chi*(p_bl_e - p_bl)/p_bl
# Interest rate on bonds - eq. 5.19
r_bl = 1/p_bl
# Expected price of bonds - eq. 5.20
p_bl_e = p_bl
# Expected capital gains - eq. 5.21
cg_e = chi*(p_bl_e - p_bl)*bl_h
# Expected regular disposable income - eq. 5.22
yd_r_e = yd_r(-1)
# Interest rate on bills - eq. 5.23
r_b = r_bar
# Price of bonds - eq. 5.24
p_bl = p_bl_bar

Let’s have a look at the graph of model LP

lp<-sfc.model(fileName="Models/LP.txt",modelName="Model Liquidity Preference")
plot.dag(lp)

The adding up of the portfolio behavior

The vertical conditions state that the sum of the exogenous components sums to 1, i.e. you cannot allocate more than the wealth and that the effect of the various return rates and disposable income are subtracting the holding of one asset to allow the increase of the holding of another.

\[\begin{align} \lambda_{10} + \lambda_{20} + \lambda_{30} &= 1 \tag{ADUP.1}\\ \lambda_{11} + \lambda_{21} + \lambda_{31} &= 0 \tag{ADUP.2}\\ \lambda_{12} + \lambda_{22} + \lambda_{32} &= 0 \tag{ADUP.3}\\ \lambda_{13} + \lambda_{23} + \lambda_{33} &= 0 \tag{ADUP.4}\\ \lambda_{14} + \lambda_{24} + \lambda_{34} &= 0 \tag{ADUP.5} \end{align}\]

The horizontal conditions, added by Godley (1996) state that the impact of an increase in the ‘own rate’ of an asset should be equal to the impact of a fall in all the other rates.

\[\begin{align} \lambda_{11} &= - \left( \lambda_{12} + \lambda_{13}\right) \tag{ADUP.6}\\ \lambda_{22} &= - \left( \lambda_{21} + \lambda_{23}\right) \tag{ADUP.7}\\ \lambda_{33} &= - \left( \lambda_{31} + \lambda_{32}\right) \tag{ADUP.8} \end{align}\]

Other authors propose the symmetry conditions that state that an increase in the return rate of one asset A will have the same impact on the holding of asset B as the increase in the return rate of asset B will have on the holding of asset A

\[\begin{align} \lambda_{12} &= \lambda_{21} \tag{ADUP.9}\\ \lambda_{13} &= \lambda_{31} \tag{ADUP.10}\\ \lambda_{23} &= \lambda_{32} \tag{ADUP.11} \end{align}\]

Note that vertical + symmetry imply automatically horizontal but that vertical + horizontal does not necessarily imply symmetry.

Higher interest rates

The first scenario run looks at the impact of interest rates (short- and long-term) on real demand. We assume that the government increases short-term rates from 3% to 4% and the long-term rate from 5% to 6.66%. We assume that this increase is completely unexpected and that the Treasury is believed in that no other change will occur.

lp<-sfc.addScenario(model=lp,vars=list(c("r_bar","p_bl_bar")),values=list(c(0.04,15)),
                            inits=1960,ends=2000)
datalp<-simulate(lp)

This replicates plot figure 5.2. p. 152

time2=c(1950:2000)
plot(time2, datalp$scenario_1[as.character(time2),"v"]/
            datalp$scenario[as.character(time2),"yd_r"],
         type="l",xlab="",ylab="",lty=1,lwd=2)
legend("bottomright",legend=c("Wealth to income ratio"),lty=c(1),lwd=2,bty="n")
grid()

We see that the increase in long-term interest rate leads to capital losses on bonds, leading to a decrease in wealth to income ratio. However because interest rates have increase, households have now a higher disposable income leading to a replenishing of their wealth. Note that as in model PC, the steady state GDP is given by

\(Y^\star=\frac{G+\left(r_B\cdot B^{\star}_h+BL_S^\star\right) \cdot (1-\theta)}{\theta} = \frac{G_{NT}}{\theta}\)

This replicates plot figure 5.3. p. 152

time2=c(1950:2000)
matplot(time2, datalp$scenario_1[as.character(time2),c("yd_r","cons")],
         type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("topleft",legend=c("Disposable income","Consumption"),lty=c(1,2),lwd=2,bty="n")
grid()

This replicates plot figure 5.4. p. 153

time2=c(1950:2000)
matplot(time2, cbind(datalp$scenario_1[as.character(time2),c("b_h")],
            datalp$scenario_1[as.character(time2),c("p_bl")]*
            datalp$scenario_1[as.character(time2),c("bl_h")])/
                datalp$scenario_1[as.character(time2),c("v")],
         type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("center",legend=c("Bonds to wealth ratio","Bills to wealth ratio"),
        lty=c(2,1),lwd=2,bty="n")
grid()

Introducing household liquidity

In order to introduce liquidity preferences, we need to modify the model. The first modification is to change the expectation on bonds price to include (i) a learning element and (ii) a stochastic term. Because expectation are going to play a role, we also need to modify the price equation for bonds such that it now reflects a corridor approach to bond pricing. The assumption is thus that the Treasury will try to pin down the price of bonds but only to a certain extent. If the demand is too large (or too low) the treasury will let prices (and interest) float. We assume that the government aims at a certain debt structure in terms of maturity.

#Adding the equations
lp2<-sfc.addEqus(lp,list(
    list(var="z1",equ = "TP>top"),
    list(var="z2",equ = "TP<bot"),
    list(var="TP",equ = "(bl_h(-1)*p_bl(-1))/(bl_h(-1)*p_bl(-1)+b_h(-1))")))
#Adding the parameters
lp2<-sfc.addVars(lp2,list(
    list(var="betae",init=0.5),
    list(var="beta",init=0.01),
    list(var="add",init=0),
    list(var="top",init=0.495),
    list(var="bot",init=0.505)
))
#Adding initial values to the endogenous variables that have lags
lp2<-sfc.editEnd(lp2,var="p_bl_e", init=20)
lp2<-sfc.editEnd(lp2,var="p_bl", init=20)
#Modifying the existing equations
lp2<-sfc.editEqus(lp2,list(
    list(var="p_bl", eq="(1+z1*beta-z2*beta)*p_bl(-1)"),
    list(var="p_bl_e",equ = "p_bl_e(-1)-betae*(p_bl_e(-1) - p_bl) + add")))
#Removing all existing scenarios (coming from LP)
lp2$scenarios<-NULL
#Checking that the model is complete
lp2<-sfc.check(lp2,fill=F)
#Adding the scenario
lp2<-sfc.addScenario(model=lp2,vars=list(c("r_bar")),values=list(c(0.035)),inits=1960,
                                            ends=2000)
lp2<-sfc.addScenario(model=lp2,vars=list(c("add")),values=list(c(-3)),inits=c(1955),
                                            ends=1955)

Let’s have a look at the graph of model LP2

plot.dag(lp2)

We are now ready to simulate the model

datalp2<-simulate(lp2)

This replicates plot figure 5.5. p. 156

time2=c(1950:2000)
plot(time2,datalp2$scenario_1[as.character(time2),"r_b"],
     lty=2, lwd=2,type="l",ylab="",xlab="",ylim=c(0.029,0.036))
par(new=T)
plot(time2, datalp2$scenario_1[as.character(time2),"r_bl"],
         type="l",xlab="",axes=F,ylab="",lty=1,lwd=2,ylim=c(0.0495,0.052))
axis(4,pretty(c(0.0495,0.052)))
legend("bottomright",legend=c("Long term interest rate",
        "Short term interest rate"),lty=c(1,2),lwd=2,bty="n")
grid()

This replicates plot figure 5.6. p. 156

time2=c(1950:2000)
matplot(time2, datalp2$scenario_1[as.character(time2),c("TP","bot","top")],
         type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend(x=1950,y=0.502,legend=c("Share of bonds in government debt held by households"),
             lty=c(2),lwd=2,bty="n")
grid()

Let’s add a second scenario by which there is an expected fall in the price of long-term bonds.

lp2<-sfc.addScenario(model=lp2,vars=list(c("add")),values=list(c(-3)),inits=c(1955),
                                            ends=1955)
datalp2<-simulate(lp2)

This replicates plot figure 5.7. p. 157

time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("r_bl")],
         type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend(x=1950,y=0.502,legend=c("Share of bonds in government debt held by households"),
             lty=c(2),lwd=2,bty="n")
grid()

This replicates plot figure 5.8. p. 158

time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("p_bl_e","p_bl")],
         type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("center",legend=c("Expected price of bonds","Actual price of bonds"),
             lty=c(1,2),lwd=2,bty="n")
grid()

This replicates plot figure 5.9. p. 158

time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("TP","bot","top")],
         type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend("center",legend=c("Share of bonds in government debt held by households"),
             lty=c(2),lwd=2,bty="n")
grid()

Modelling the firm: price, investment and profits

The role of banks: credit creation and investment

Functional income distribution

Consumption function

Steady state of BMW

#Supply of consumption goods - eq. 7.1
 c_s = c_d
#Supply of investment goods - eq. 7.2
 i_s = i_d
#Supply of labour - eq. 7.3
 n_s = n_d
#Transactions of the firms
#GDP - eq. 7.5
 y = c_s + i_s
#Wage bill - eq. 7.6
 wb_d = y - r_l(-1)*l_d(-1) - af
#Depreciation allowances - eq. 7.7
 af = delta*k(-1)
#Demand for bank loans - eq. 7.8
 l_d = l_d(-1) + i_d - af
#Transactions of households
#Disposable income - eq. 7.9
 yd = wb_s + r_m(-1)*m_h(-1)
#Bank deposits held by households - eq. 7.10
 m_h = m_h(-1) + yd - c_d
#The wage bill
#"Supply" of wages - eq. 7.13
 wb_s = w*n_s
#Labour demand - eq. 7.14
 n_d = y/pr
#Wage rate - eq. 7.15
 w = wb_d/n_d
#Household behaviour
#Demand for consumption goods - eq. 7.16
 c_d = alpha0 + alpha1*yd + alpha2*m_h(-1)
#The investment behaviour
#Accumulation of capital - eq. 7.17
 k = k(-1) + i_d - da
#Depreciation allowances - eq. 7.18
 da = delta*k(-1)
#Capital stock target - eq. 7.19
 k_t = kappa*y(-1)
# below is an additional equation I've defined to get output to capital ratio as in figure 7.4, the model will run even without equation
OCR = y/k(-1)  
#Demand for investment goods - eq. 7.20
 i_d = gamma*(k_t - k(-1)) + da

Out of equilibrium values

BMW<-sfc.model("Models/BMW.txt")
## Warning in sfc.check(model, fill = fill): The following variables have lags
## but no initial values: - r_l - r_m
plot.dag(BMW)

Stability

More on profits and inventories

Components Firms Current account Firms Capital account
Sales \(+S\)
Change in the value of inventories +\(\Delta IN\) -\(\Delta IN\)
Wages \(-WB\)
Interest on loans \(-r_{l-1} L_{-1}\)
Entrepreneurial profits \(-F\)
Change in loans \(+\Delta L\)
Sum 0 0

Calibration/Estimation

We will now briefly discuss the different ways to calibrate/estimate a model and then try the various solutions to a (extremely) simple model: PCEX.

Steady/stationary state calibration

Assuming you want to start with a model at the steady/stationary state, one way of calibrating these model is by doing “forward deduction”:

  1. Find in the literature (SFC or other) the values for the parameters you have in your model
  2. Run the model, hoping it will converge towards the steady state
  3. Fiddle as needed with some parameters such that the steady state meets your requirements

Steady/stationary state calibration

Assuming you want to start with a model at the steady/stationary state, one way of calibrating these model is by doing “backward induction”:

  1. First you need to change your model equation so that it represent a steady state model
    1. All lags are removed
    2. All expectation are equal to the corresponding steady state values
    3. Some equation have to be changed completely: i.e. assuming a change in wealth function (\(V=V_{-1}+YD-C\)), it becomes (\(C=YD\)) but then \(V\) is undefined. This implies that the consumption function (\(C = \alpha_1 YD^e + \alpha2 V_{-1}\)) has to change as well so as to show the relationship between wealth and income (\(V=\frac{1-\alpha_1}{\alpha_2}YD\))
  2. Endogenise the parameters and exogenise some stock-flow norms or accepted values for specific indicators
    1. Endogenising the parameters imply that instead of having the wealth to income function (\(V=\frac{1-\alpha_1}{\alpha_2}YD\)), we now specify \(\alpha_2\) as being determined by this equation: \(\alpha_2 = (1-\alpha1)\frac{YD}{V}\).
    2. Exogenising the stock-flow norms or accepted indicators implies fixing debt-to-gdp ratios or unemployment rates to realistic values
  3. Once you have enough equations for the number of parameters, and enough constraints (i.e. fixed values for stock-flow norms), the model should be determined and it is possible to solve the system and obtain both the parameters and initial values for the specific steady state. Note that you can use the PK-SFC package for that.

Out-of-Steady State calibration/estimation: massaged data

one of the problem with SFC models is that they explicitly constrain the data into a specific data structure. I.e. the balance sheets and other transaction flow matrices are sparse while it is not the case in reality. This means that you will have to massage the data so that it meets your requirement. There is a lot of subjectivity happening there. A pseudo-algorithm for this part would look like this.

  1. Individuate the time series you want to be exactly like the real world (e.g. GDP, unemployment, Household wealth)
  2. Individuate the time series you are ok to play the buffer: these will have to be constructed from the accounting structure of the model. E.g. if you have households having two type of assets: money and bonds and you want households wealth and bond holding to be based on data, then money will have to be constructed as wealth-bonds.
  3. Once you have a stock-flow consistent data set: run regressions to compute the parameters
    1. I personally use OLS or dynamic OLS equation per equation
    2. You could use 2-stage LS or minimum likelihood to estimate the system as a whole
    3. Gennaro Zezza and others use co-integration methods
  4. Initial conditions are determined by the data you have.
  5. You need to say what you will do with residuals (coming out of your regressions): set them a 0, have them fixed at a non-zero value, have them decreasing through time?

Pure estimation

If you want to do a pure empirical model, then you have to start from the data and construct the model around it. You will need to look at how disaggregated the data is and what you can afford to do. A pseudo algorithm would look like this:

  1. Aggregate/Disaggregate existing data sets such as it meets with the model structure you think is relevant or modify the model structure you want to have such that it fits with the data you have
  2. Once you are happy with the data and the model structure, run regressions to compute the parameters
    1. I personally use OLS or dynamic OLS equation per equation
    2. You could use 2-stage LS or minimum likelihood to estimate the system as a whole
    3. Gennaro Zezza and others use co-integration methods
  3. Initial conditions are determined by the data you have.
  4. You need to say what you will do with residuals (coming out of your regressions): set them a 0, have them fixed at a non-zero value, have them decreasing through time?

Practical example

We will try to apply the middle two techniques (Backward induction and massaged data) to model PCEX to the UK economy.

Backward induction

library(PKSFC)

# Loading the Steady State version of the model: you need to write down that file. 
# Bear in mind that you dont need to have more than two periods in the steady state
modelSS<-sfc.model("Models/PCEXSS.txt")

plot.dag(modelSS)

# Solving the steady state
dataSS<-simulate(modelSS)
View(dataSS$baseline)
#Getting the value for all parameters and steady state values
calib<-as.data.frame(t(dataSS$baseline[2,]))

#What follows is a bit tricky but not too complicated to understand. What we will do is 
# to create a temporary file called calibration.txt based on the values obtained via the
# steady state model. This file will then be merged into the PCEX_forCalibration.txt file 
# so as to dynamically create a complete model file (i.e. equations + calibration) which 
# can then be simulated

#Creation of the calibration.txt file
sink("calibration.txt")
#######  PARAMETERS
cat(paste("alpha1=",calib$alpha1,"\n"))
cat(paste("alpha2=",calib$alpha2,"\n"))
cat(paste("lambda0=",calib$lambda0,"\n"))
cat(paste("lambda1=",calib$lambda1,"\n"))
cat(paste("lambda2=",calib$lambda2,"\n"))
cat(paste("theta=",calib$theta,"\n"))
#######  EXOGENOUS
cat(paste("g=",calib$g,"\n"))
cat(paste("r_bar=",calib$r_bar,"\n"))
#######  Starting values for stocks
cat(paste("b_cb=",calib$b_cb,"\n"))
cat(paste("b_h=",calib$b_h,"\n"))
cat(paste("b_s=",calib$b_s,"\n"))
cat(paste("h_h=",calib$h_h,"\n"))
cat(paste("h_s=",calib$h_s,"\n"))
cat(paste("r=",calib$r,"\n"))
cat(paste("v=",calib$v,"\n"))
cat(paste("yd=",caliby$d,"\n"))
sink()

# Now we will read the calibration file
calibrationLines<-readLines("calibration.txt")

#Next step is to read the model without calibration
modelLines<-readLines("PCEX_forCalibration")

#We need to determine where to add the calibration lines, this is done where 
# "{CALIBRATION}" has been inserted in the PCEX_forCalibration.txt file
indexCalib<-grep("{CALIBRATION}",modelLines)

#Now we can merge the two files into a complete model
totModel<-c(modelLines[1:(indexCalib-1)],calibrationLines,
                        modelLines[(indexCalib+1):length(modelLines)])

#Write down the model into a txt file
writeLines(totModel,"Models/PCEX_Auto.sfc")

#And load the model. The model is now ready to be simulated
WholeModel<-sfc.model("Models/PCEX_Auto.sfc")

Massaged data

library(PKSFC)

#1. fetch all the relevant data sets for your case, using pdfetch
#getting the data in a raw format
NF_TR_raw = pdfetch_EUROSTAT("nasa_10_nf_tr", UNIT="CP_MNAC", NA_ITEM=c("..."),
                                            SECTOR=c("..."), DIRECT = "...")
#Transforming the data into a data.frame (easy to access and to view)
NF_TR<-as.data.frame(NF_TR_raw)

#2. Construct the various massaged time series
NF_TR$MASSAGED<-NF_TR[,"..."]+NF_TR[,"..."]-NF_TR[,"..."]

#3. Run the regression
model.cons<-lm(CONS[-1]~YD[-19]+V[-19]-1,data=NF_TR)
calib$alpha1<-coef(model.cons)[1]
calib$alpha2<-coef(model.cons)[2]

...

#What follows is a bit tricky but not too complicated to understand. What we will do is 
# to create a temporary file called calibration.txt based on the values obtained via the
# steady state model. This file will then be merged into the PCEX_forCalibration.txt file 
# so as to dynamically create a complete model file (i.e. equations + calibration) which 
# can then be simulated

#Creation of the calibration.txt file
sink("calibration.txt")
#######  PARAMETERS
cat(paste("alpha1=",calib$alpha1,"\n"))
cat(paste("alpha2=",calib$alpha2,"\n"))
cat(paste("lambda0=",calib$lambda0,"\n"))
cat(paste("lambda1=",calib$lambda1,"\n"))
cat(paste("lambda2=",calib$lambda2,"\n"))
cat(paste("theta=",calib$theta,"\n"))
#######  EXOGENOUS
cat(paste("g=",calib$g,"\n"))
cat(paste("r_bar=",calib$r_bar,"\n"))
#######  Starting values for stocks
cat(paste("b_cb=",calib$b_cb,"\n"))
cat(paste("b_h=",calib$b_h,"\n"))
cat(paste("b_s=",calib$b_s,"\n"))
cat(paste("h_h=",calib$h_h,"\n"))
cat(paste("h_s=",calib$h_s,"\n"))
cat(paste("r=",calib$r,"\n"))
cat(paste("v=",calib$v,"\n"))
cat(paste("yd=",caliby$d,"\n"))
sink()

# Now we will read the calibration file
calibrationLines<-readLines("calibration.txt")

#Next step is to read the model without calibration
modelLines<-readLines("PCEX_forCalibration")

#We need to determine where to add the calibration lines, this is done where
# "{CALIBRATION}" has been inserted in the PCEX_forCalibration.txt file
indexCalib<-grep("{CALIBRATION}",modelLines)

#Now we can merge the two files into a complete model
totModel<-c(modelLines[1:(indexCalib-1)],calibrationLines,
                    modelLines[(indexCalib+1):length(modelLines)])

#Write down the model into a txt file
writeLines(totModel,"Models/PCEX_Auto.sfc")

#And load the model. The model is now ready to be simulated
WholeModel<-sfc.model("Models/PCEX_Auto.sfc")

References

Foley, D. 1975. “On Two Specifications of Asset Equilibrium in Macroeconomic Models.” Journal of Political Economy 83 (2).

Godley, Wynne. 1999. “Seven Unsustainable Processes: Medium-Term Prospects and Policies for the United States and the World.” Strategic Analysis. The Levy Economic Institute of Bard College.